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ABSTRACT 

High order essentially non-oscillatory (ENO) schemes, originally designed for compress- 
ible flow and in general for hyperbolic conservation laws, are applied to incompressible 
Euler and Navier-Stokes equations with periodic boundary conditions. The projection to 
divergence-free velocity fields is achieved by fourth order central differences through Fast 
Fourier Transforms (FFT) and a mild high-order filtering. The objective of this work is to 
assess the resolution of ENO schemes for large scale features of the flow when a coarse grid is 
used and small scale features of the flow, such as shears and roll-ups, are not fully resolved. 
It is found that high-order ENO schemes remain stable under such situations and quantities 
related to large-scale features, such as the total circulation around the roll-up region, are 
adequately resolved. 
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1 Introduction 


In this paper we consider numerically solving the incompressible Navier-Stokes or Euler 
equations 


lit H" HHx VUy — pi^xx ^yy) Px 

V t + UV X + VVy = fl(V XX + Vyy) - Py (l.l) 

1l X “I” Vy — 0 

or their equivalent conservative form 

U t + ( U 2 )x + (uv)y = p(u xx + Uyy) - p x 

Vt T jj “I - (u — p[v xx T Uyy) py (1.2) 

U X + Vy = 0 

where (u, v) is the velocity vector, p is the pressure, p > 0 for the Navier-Stokes equations and 
p = 0 for the Euler equations. The numerical methods we use are the high-order essentially 
non-oscillatory (ENO) schemes, originally designed for compressible flow and in general for 
hyperbolic conservation laws [8], [12]. The equation is defined on the box [0,2 tt] x [0,2tt] 
with periodic boundary conditions in both directions. We choose two space dimensions for 
easy presentation, although our method is also applicable for three space dimensions. 

In some sense equations (1.1) are easier to solve numerically than their compressible 
counter-parts because the latter have solutions containing possible discontinuities (for ex- 
ample shocks and contact discontinuities). However, the solution to (1.1), even if for most 
cases smooth mathematically, may evolve rather rapidly with time t and may easily become 
too complicated to be fully resolved on a feasible grid. Traditional linearly stable schemes, 
such as spectral methods and high-order central difference methods, are suitable for the 
cases where the solution can be fully resolved, but typically produce signs of instability such 
as oscillations when small scale features of the flow, such as shears and roll-ups, cannot be 
adequately resolved on the computational grid (see, for example, the last contour in Fig. 1). 
Although in principle one can always overcome this difficulty by refining the grid, today’s 
computer capacity seriously restricts the largest possible grid size. 

In the last few years there is a lot of activity in designing high-order, nonlinearly stable 
“shock capturing” schemes for compressible flow and in general for hyperbolic conservation 
laws. See, for example, [2], [7], [8], [12], and the references listed therein. The philosophy of 
such schemes is to give up fully resolving rapid transition regions or shocks, just to “capture” 
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them in a stable and somehow globally correct fashion (e.g., with correct shock speed), but at 
the same time to require a high resolution for the smooth part of the flow. The success of such 
an approach for the compressible flow is documented by many examples in the literature. 
One example is the one and two dimensional shock interaction with vorticity or entropy 
waves [12], [13]. The shock is captured sharply and certain key quantities related to the 
interaction between the shock and the smooth part of the flow, such as the amplification and 
generation factors when a wave passes through a shock, are well resolved. Another example is 
the homogeneous turbulence for compressible Navier-Stokes equations studied in [13]. In one 
of the test cases, the spectral method can resolve all the scales using a 256 2 grid, while third 
order ENO with just 64 2 points can adequately resolve certain interesting quantities such as 
the average fluctuation Mach number and maximum divergence, although it cannot resolve 
local quantities achieved inside the rapid transition region such as the minimum divergence. 
The conclusion seems to be that, when fully resolving the flow is either impossible or too 
costly, a “capturing” scheme such as ENO can be used on a coarse grid to obtain at least 
some partial information about the flow. 

In this paper we perform a similar numerical study for the incompressible equation (1.1), 
to see what one can get by using the high-order ENO schemes on a coarse grid, without 
fully resolving the flow. We choose doubly periodic shear layers, used in [1], as our test case. 
A spectral method with 512 2 points seems able to fully resolve the flow up to t = 8, but 
begins to show signs of underresolution (wriggles in vorticity) thereafter. This indicates the 
tremendous requirement upon computation resources if one tries to resolve everything in the 
flow. We then use the third order ENO method (which is fourth order in the Lj sense) and 
coarse grids (64 2 and 128 2 points), to assess its resolution. We find that the scheme remains 
stable for coarse meshes and certain quantities related to the smooth part of the flow, such 
as the total circulation around the roll-up, are adequately resolved by ENO methods. 

A pioneer work in applying shock capturing compressible flow techniques to incompress- 
ible flow is by Bell, Colella and Glaz [1], in which they considered a second order Godunov 
type discretization, investigated the projection into divergence-free velocity fields for general 
boundary conditions, and discussed accuracy of time discretizations. Since our objective in 
this paper is to assess the resolution of the ENO method for (1-1), we choose a periodic 
boundary condition to simplify the projection. General boundary conditions would have to 
be handled either by more complicated projection [l] or by artificial compressibility methods 
[4]. We are currently investigating ENO schemes for such cases. 
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2 The ENO Method 


We solve (1.2) in its equivalent projection form 


u 

v 


= P 


u 

uv 



+ 


u 

V 


yyj 


where P is the Hodge projection into divergence-free fields, i.e., if 


u 

v 


( 2 . 1 ) 


| , then 


u x + v y = 0 and v y — u x = v y — u x . See, e.g., [1]. For the current periodic case the additional 
condition to obtain a unique projection P is that the mean values of u and v are preserved, 
be., Jo* Jo 2 * “(*» y)dxdy = JjJ* «(*, y)dxdy and / 0 2x / 0 2x u(x, t/)dxdy = / 0 2x / 0 2x v(x, y)dxdy. 

We use N x and N y (even numbers) equally spaced grid points in x and j/, respectively. 
The grid sizes are denoted by Ax = ^ and Ay = and the grid points are denoted by 
X{ = iAx and y 3 = jAy. The approximated numerical values of u and v at the grid point 
(xi,yj) are denoted by U{j and Vij. 

We first describe the numerical implementation of the projection P. In the periodic case 
this is easily achieved in the Fourier space. We first expand u and v using Fourier collocation: 


Ujl 

2 


Ejl 

2 


U N 


(*,y) = Z Z 


Ukie 


I(kx+ly) 


UlL 

2 


N x 

2 


v N (x,y)= Z Z 


v kl e 


I (kx+ly) 


(2.2) 


l=-3“= 


Nx 

2 




where / = y/—l, u k i and v k i are the Fourier collocation coefficients which can be computed 
from the point values Uij and using either FFT or matrix- vector multiplications. The de- 
tail can be found in, e.g., [3]. Derivatives, either by spectral method or by central differences, 
involve only multiplications by factors d k or df in (2.2) because e / ( <:x+,y ) are eigenfunctions 
of such derivative operators. For example, 


d x k = Ik , df - II 


for spectral derivatives; 


2/sin(lf0 


jx _ 2/ sin( 2 X ) _ — v 2 

dk ~ Z ’ d ‘ ~ aT“ 


(2.3) 


(2.4) 


for the second order central differences which, when used twice, will produce the second 
order central difference approximation f or w xx , and 


a k ~ 


21 yj (1 — cos(A:Ax))(7 — cos(A:As)) ^ 2/^/(l — cos(/Ay))(7 — cos (/Ay)) 

~Kx ’ d ‘ = A y 


(2.5) 
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for the fourth order central differences which, when used twice, will produce the fourth order 
central difference approximation 1( ^" +] - W| ~ 1 ■ g* 2 + for w XT . High order filters, such 

as the exponential filter [10], [9]: 


,-“(7fc) 2p 


y -«(-sr) 2p 
erf — e N v ’ 


- o » (2.6) 

where 2 p is the order of the filter and a is chosen so that e~ a is machine zero, can be used 
to enhance the stability while keeping at least 2p-th order of accuracy. This is especially 
helpful when the projection P is used for the under- resolved coarse grid with ENO methods. 
We use the fourth order projection (2.5) and the filter (2.6) with 2p = 8 in our calculations. 
This will guarantee third order accuracy (fourth order in L x ) of the ENO scheme. We will 
denote this combination (the fourth order projection plus the eighth order filtering) by P 4 . 

To be precise, if^*?^=P 4 ^^^ and u kX and v m are Fourier collocation coefficients of u 
and v, then the Fourier collocation coefficients of u and v are given by 


u 


Wt 


yd^jc^u - d%v) 


(d x k y + K ) 2 3 ’ u ak<7 ‘ {diy + m 2 

where cr x k and erf are defined by (2.6) with 2 p = 8, and d x k and df are defined by (2.5). 

Next we shall describe the ENO scheme for (2.1). Since (2.1) is equivalent to the non- 
conservative form (1.1), it is natural to implement upwinding by the signs of u and u, and 
to implement ENO equation by equation. The r-th order ENO approximation of, e.g., ( u 2 ) x 
is thus summarized as follows: 


-d k (tfu- d%v) 


(2.7) 


1. Take f(x) = u 2 (x,y) with y fixed. We start with the point values /,• = /(s,); 

2. Define a function h(x) satisfying f(x) = f x _dl h(Qd(, and its primitive H(x) = 
f x h(£)d(. For each j + i, a (r + l)-th order polynomial Q j+x / 2 (x) is constructed 
interpolating H(x) at points near Xj + 1 / 2 . It is remarked in [12] that the Newton 
divided differences of H(x) are constant multiples of those of f(x) of one order lower. 
Since we do not need the zeroth order difference, the approximation Q J+ 1/2 of H(x) 
can be accomplished using the information of fi only, without explicitly constructing 
H(x) or its differences; 

3. The stencil of the polynomial Q(x) is determined adaptively by upwinding and smooth- 
ness of f(x). It starts with either Xj or x j+ x according to whether u > 0 or u < 0, then 
one point to the left or right is added to the stencil each time by comparing the two 
relevant divided differences and picking the smaller one in magnitude. 
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4. The derivative f'(x ) is finally approximated by the conservative difference -^{fj+ 1/2 — 
fj- 1 / 2 ) where the numerical flux is computed by Jj+ 1/2 = d x Q( x ) W=^ J+ i/2 • 

5. In the actual coding of the algorithm, undivided differences should be used instead of 
the divided differences to reduce round-off errors. There are also ways to make the 
procedure more economical on a vector computer. The details can be found in [13]. 

The approximation to ( uv) x is the same as above with f(x) = u(x,y)v(x,y), and that 
for (uv) y and for ( v 2 ) y are similar, with upwinding based on v. 

There are two ways to handle the second derivative terms for the Navier-Stokes equations. 
One can absorb them into the convection part and treat them using ENO. For example, 
f(x) = u 2 (x,y) can be replaced by f(x) = u 2 (x,y ) - pu(x, y) x , where u(x,y) x itself can be 
obtained using either ENO or central difference of a suitable order. The remaining procedure 
for computing /(x) r would be the same as described above. Another simpler possibility is 
just to use standard central differences (of suitable order) to compute the double derivative 
terms. Our experience with compressible flow is that there is little difference between the 
two approaches, especially when the viscosity n is small. 

In the above we have described the discretization for the spatial derivatives 



y = yj 


Next we use a third order TVD (total variation diminishing) Runge-Kutta method [11] to 
discretize the resulting ODE: 

(“) ( = p ,i, (2.9) 

obtaining: 



Notice that we have used the property P 4 0 P 4 = P 4 in obtaining the discretization (2.10) 
from (2.9). 
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This explicit time discretization is expected to be nonlinearly stable under the CFL 
condition 


At 


max 

hi 


Kl + M ) + 2 y 


i i 

+ 


Ax 2 Ay' 2 1 


< 1 


( 2 . 11 ) 


Ax Ay 

(see [1 1]). For small y (which is the case we are interested in) this is not a serious restriction 
on At. 


3 Numerical Results 

We present two numerical examples in this section. 

Example 1 : This example is used to check the third order accuracy of our ENO scheme 
for smooth solutions. We first take the initial condition as 

u(x, y, 0) = — cos(x) sin(y), u(x, y, 0) = sin(x) cos(y) (3-1) 

which was used in [5]. The exact solution for this case is known: 

u(x, y, t) — — cos(x) sin(y)e _2M< , v(x,y,t) = sin(x)cos(y)e~ 2fit (3-2) 

We take Ax = Ay = with N = 32,64, 128 and 256. The solution is computed up 
to t = 2 and the L 2 error and numerical order of accuracy are listed in Table 3.1. For the 
y = 0.05 case, we list results both with fourth order central approximation to the double 
derivative terms (central) and with ENO to handle the double derivative terms by absorbing 
them into the convection part (ENO). We can clearly observe fully third order accuracy 
(actually better in many cases because the spatial ENO is fourth order in the L\ sense) in 
this table. 


Table 3.1: Accuracy of ENO Schemes for (3.1) 


N 

y = 0 

y — 0.05, central 

y = 0.05, ENO 

Z /2 error 

order 

Z .2 error 

order 

L 2 error 

order 

32 

9.10(-4) 


5.28(-4) 


4.87(-4) 


64 

5.73(-5) 

3.99 

3.20(-5) 

4.04 

3.09(-5) 

3.98 

128 

3.62(-6) 

3.98 

1.93(-6) 

4.05 

1.89(-6) 

4.03 

256' 

2. 28 (-7) 

3.99 

1.18(-7) 

4.03 

1.16(-7) 

4.03 


We then take the initial condition as 
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u(x,y,0) — — sin 2 (x) sin(2y), *>(x,y,0) = sin(2x) sin 2 (y) (3.3) 

which was used in [1]. We again take Ax = Ay — ^ with N = 16,32,64, 128 and 256 and 
compute the solution up to t = 2. However, this time the exact solution is unknown. As 
in [1], we test the accuracy by computing the L 2 difference between the solutions at the 
grid sizes 2 Ax and Ax, on the coarser grid. If w& x = w + CAx T + 0(Ax r+1 ), then this 
difference w 2 &x ~ w &x = (2 r — l)CAx T + 0(Ax r+1 ) would predict the correct asymptotic 
order of accuracy r and the error itself on the finer grid multiplied by 2 r — 1. The result is 
summarized in Table 3.2, where w = ( u,v) T and L 2 diff= ||w 2 Ax — w Axi|L 2 - 1° 4hi s table 
the estimated order and error are obtained using the remarks above. For p = 0.05, we have 
again provided results both with fourth order central approximation to the double derivative 
terms (central) and with ENO to handle the double derivative terms by absorbing them 
into the convection part (ENO). The result in Table 3.2 reconfirms better than third order 
accuracy for the third order ENO scheme. 

Table 3.2: Accuracy of ENO Schemes for (3.3) 


N 

p = 0 

p = 0.05, central 

p = 0.05, ENO 

L 2 diff 

order 

error 

L 2 diff 

order 

error 

L 2 diff 

order 

error 

32 

i-u(-i) 



3.20(-2) 



3.60(-2) 



64 

1.40(-2) 

3.02 

1.96(-3) 

2.78(-3) 

3.52 

2.66(-4) 

2.93(-3) 

3.62 

2.60(-4) 

128 

1.46(-3) 

3.26 

1.69(-4) 

1.81 (-4) 

3.94 

1.26(-5) 

1.80(-4) 

4.02 

1.18(-5) 

256 

1-11 (-4) 

3.77 

8.78(-6) 

1.09(-5) 

4.06 

6.91(-7) 

1.10(-5) 

4.04 

7.15(-7) 


Example 2: This is our test example to study resolution of ENO schemes when the grid 
is coarse. It is a double shear layer taken from [1]: 


u(x,y,0) = <Ssin(x) 


(3.4) 


ll(r „ = / tanh((y - tt/2 )/p) y< n 

' ( tanh((37r/2 — y)/p) y>n 

where we take p — 7r/15 and 8 = 0.05. The Euler equations (p = 0) are used for this example. 
The solution quickly develops into roll-ups with smaller and smaller scales, so on any fixed 
grid the full resolution is lost eventually. For example, the expensive run we performed using 
512 2 points for the spectral collocation code (with a 18-th order filter (2.4)) seems able to 
resolve the solution fully up to t = 8, then begins to lose resolution as indicated by the 
wriggles in the vorticity contour at t = 10 (the last contour in Fig. 1). On the other hand, 
the ENO runs with 64 2 and 128 2 points produces smooth, stable results (Fig. 2 and 3). 
Apparently with these coarse grids the full structure of the roll-up is not resolved. However, 
when we compute the total circulation 
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(3.5) 


cq— uj(x,y)dxdy = / udx + vdy 
Jq Jdu 

around the roll-up by taking H = [f, x [0,2 tt] and using the rectangular rule (which is 
infinite order accurate for the periodic case) on the line integrals at the right-hand-side of 
(3.5), we can see that this number is resolved much better than the roll-up itself (see Table 
3.3). 


Table 3.3: Resolution of the Total Circulation 


t 

2 

4 

6 

8 

10 

ENO 64 2 

0.87300 

3.07100 

7.16889 

9.88063 

10.90122 

ENO 128 2 

0.87452 

2.97810 

7.30999 

10.34414 

11.79418 

spectral 51 2 2 

0.87433 

2.98029 

7.28308 

10.46212 

11.85875 


4 Concluding Remarks 

High order ENO schemes can be applied to the incompressible Euler and Navier-Stokes 
equations to obtain stable, under-resolved results on a coarse grid. Some global quantities 
such as the circulation around the roll-up region can be faithfully resolved. 
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